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We investigate the properties of the Ising- Glauber model on a periodic cubic lattice of linear 
dimension L after a quench to zero temperature. The resulting evolution is extremely slow, with 
long periods of wandering on constant energy plateaux, punctuated by occasional energy-decreasing 
spin- flip events. The characteristic time scale r for this relaxation grows exponentially with the 
system size; we provide a heuristic and numerical evidence that r ^ exp(L^). For all but the 
smallest-size systems, the long-time state is almost never static. Instead the system contains a 
small number of "blinker" spins that continue to flip forever with no energy cost. Thus the system 
wanders ad infinitum on a connected set of equal-energy blinker states. These states are composed 
of two topologically complex interwoven domains of opposite phases. The average genus ql oi the 
domains scales as L^, with 7 1.7; thus domains typically have many holes, leading to a "plumber's 
nightmare" geometry. 

PACS numbers: 64.60.My, 75.40.Gb, 05.50.+q, 05.40.-a 



I. INTRODUCTION 

We study the evolution of the homogeneous Ising fer- 
romagnet on a periodic cubic lattice in which the spins 
are endowed with zero- temperature Glauber dynamics. 
Starting from an initial state of zero magnetization, cor- 
responding to a supercritical temperature, conventional 
wisdom states that the spins organize into a coarsening 
domain mosaic whose characteristic length scale grows as 
t^/^ [ll-Q- This coarsening continues as long as the typ- 
ical domain size is less than the linear dimension of the 
system L. At longer times, it is natural to anticipate that 
one of the ground states, with magnetization m = +1 or 
m = — 1, should ultimately be reached. 

For the one-dimensional system, this expectation is ful- 
filled — the ground state is always reached, independent 
of the initial condition. As one might expect by the dif- 
fusive dynamics of the spin domain interfaces, the char- 
acteristic time to reach this final state scales as L^. In 
two dimensions, the ground state is no longer the only 
asymptotic outcome. A system that starts at zero mag- 
netization may also get stuck in an infinitely long-lived 
metastable state that consists of straight single-phase 
stripes [s"-^ . The probability to reach such a stripe state 
was found numerically to be close to |. Recently, a the- 
oretical argument was given [9] that relates the stripe 
state probability to certain exactly-calculated percola- 
tion crossing probabilities. The resulting probabilities to 
reach the stripe state are ^ ~ 2^ f| = 0.3558 ... for free 
boundary conditions and 0.3390 ... for periodic boundary 
conditions, in agreement with simulation data @, @. 

For regular lattices in three dimensions or greater, little 
is known about the evolution and long-time state of this 
kinetic Ising ferromagnet. If the initial magnetization 
is non-zero, it is generally believed that this system (on 
an even-coordinated lattice in d > 2) ultimately reaches 
the ground state of the majority phase in the thermo- 
dynamic limit, no matter how small the initial magne- 
tization [10]. This result has not been proved, even in 



two dimensions, although it is plausible because of the 
connection with percolation [9]. In the d ^ 00 limit, the 
fact that the (majority) ground state is reached is intu- 
itively obvious and has been recently proved in Ref. [lOj]. 
Physically, however, the zero initial-magnetization state 
is much more important than the general case of a non- 
zero magnetization. Indeed, the usual coarsening process 
begins at a temperature that exceeds the critical tem- 
perature where the initial magnetization (of an infinite 
system) equals zero. We therefore focus on the case of 
zero initial magnetization in this paper. 

An earlier study [6*] found that an L x L x L Ising ferro- 
magnet on a cubic lattice exhibits a much more complex 
evolution than the corresponding two-dimensional sys- 
tem. In particular, the long-time state is typically topo- 
logically complex and not static. Here, we investigate the 
three-dimensional system in greater detail [11]. We offer 
new perspectives to help understand the long-time state 
of the system and present simulation results to quantify 
its unusual properties. 

While the ground state is a possible long-time outcome 
of the dynamics, it happens that a geometrically rich and 
infinitely long-lived metastable state (Fig. [T]) arises with 
overwhelming probability. This long-time state typically 
has a sponge-like geometry, with multiple interpenetrat- 
ing regions of positive and negative magnetization. The 
continuum version of this state is one in which the mean 
curvature of the interface is zero. This restriction leads to 
a veritable zoo of possi ble g eometries that have been ex- 
tensively cataloged [HI, |lj A few illustrative examples 
of a subclass of these systems — periodic zero-curvature 
continuum interfaces — are given in Fig. [2l These struc- 
tures also resemble the geometrically complex arrange- 
ments that arise in two-phase micellar systerns, and are 
popularly known as "plumber's nightmares" jl4l4l7j . 

Intriguingly, the long-time state of the three- 
dimensional system is generally not static, but rather, 
contains stochastic blinker spins that can flip repeatedly 
without any energy cost. The interfaces defined by these 
spins can therefore wander ad infinitum on a connected 
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FIG. 1. (color online) A typical example of a state with 
blinker spins for a 20^ cubic lattice with periodic boundary 
conditions. The highlighted blocks indicate blinker spins. 



pendixIEl Typically, domain interfaces have large genus 
so that the domains have many interpenetrating protru- 
sions. Finally, we investigate the time dependence of the 
survival probability (Sect. HVD]) . namely, the probability 
that the energy of the system is still decreasing up to a 
given time. Based on the insights gained from studying 
blinker states, we can understand some important fea- 
tures of this survival probability. Concluding remarks 
are given in Sect. [Vl Basic features of the evolution of 
a 2^ system, where all states can be enumerated exactly, 
as well as a few details of slightly larger systems are are 
presented in Appendix [Bl 



II. MODEL 



The Hamiltonian of the ferromagnetic Ising model is 
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FIG. 2. Triply-periodic minimal surfaces in three dimensions. 
Figures provided by K. Brakke (see also Ref. [17i]). 



but bounded set of equal-energy blinker states as illus- 
trated in Fig. [TJ Here the cubes correspond to up spins 
(with the spin at the cube center), while the down spins 
correspond to blank space. The highlighted cubes indi- 
cate spins at the convex (outer) corners of domain in- 
terfaces that can flip with no energy cost. There are 
also oppositely-oriented spins (blank spaces) adjacent to 
the apex of the concave corners that can also flip up 
with no energy cost. A blinker state wanders perpetu- 
ally on a small set of iso-energy points in state space. 
Blinker states first appear (albeit rarely) when the linear 
dimension L = 5, but essentially all configurations con- 
tain blinkers for large L. While the fraction of blinker 
spins is quite small, the fraction of the system volume 
over which blinker spins can wander is macroscopic — it 
is usually of the order of ten percent for large L. 

In Sect, ini we define the system under study. In 
Sect. we present a simple method to accelerate the 
simulations. Details about the accuracy of this accel- 
eration algorithm are given in Appendix [A] We then 
discuss the physical properties of the long-time state of 
the system in Sect. \IV\ including details about blinker 
states (Sect. IIV Ap . the asymptotic energy of the sys- 
tem (Sect. lIVBp . and topological characteristics of the 
domains (Sect. lIVCp . Additional facts about the L- 
dependence of some basic observables are given in Ap- 



where <Ji = ±1 denotes the spin at site z, the interaction 
strength is set to one, and the sum is over all nearest- 
neighbor pairs of sites. If not stated explicitly other- 
wise, we consider the cubic lattice of linear dimension 
L, with L even, and with periodic boundary conditions. 
There are two natural choices for the initial state: (i) 
initially uncorrelated and equal fractions of +1 and — 1 
spins (corresponding to an initial temperature T = oo), 
or (ii) the antiferromagnetic initial condition, in which all 
pairs of neighboring spins are oppositely oriented. The 
long-time evolutions of the system starting from these 
two initial states are similar, with only small quantita- 
tive differences in the distribution of basic physical ob- 
servables, such as the energy and magnetization in the 
long-time state. Thus we focus on the antiferromagnetic 
initial state for concreteness and for simplicity. 

The spins evolve by zero-temperature single spin-flip 
Glauber dynamics [18]. To implement this dynamics at 
zero temperature, we keep a hst of flippable spins — those 
where the energy change of the system AE would be 
zero or negative if the spin were to flip. (At zero tem- 
perature, spin- flips which would lead to an increase of 
energy, AE > 0, are forbidden [^.) By picking only 
the spins from this list, we eliminate the time that would 
be wasted in picking and simulating no n- flippable spins. 
In each update, we pick a flippable spin at random and 
flip it with probability 1 if AE < 0, or with probabil- 
ity ^ if = 0. This update corresponds to majority 
rule, as the condition AE < means that the majority 
of neighbors are antiparallel to the selected spin. 

More generally, zero-temperature single spin-flip dy- 
namics is defined by the rules. 



Flipping probability = < 



if AE < 0, 
if AE = 0, 
if AE>0, 



(2) 
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that depend on the single parameter p. For the heat-bath 
and Glauber dynamics, p = |, while for the Metropolis 
algorithm [l9| all updates proceed with the same rate 
{p = 1). While the behavior for all p > is essentially 
independent of the p = case is special, as the system 
quickly gets trapped in a jammed configuration [20|. We 
shall avoid the p = dynamics with only strictly energy- 
lowering spin fiips, as this case is relevant for a special 
class of kinetically constrained systems (see ^ for re- 
view). Since the choice for the parameter p is arbitrary 
(as long as p is strictly positive), we fix p = |, although 
this choice is essentially a matter of habit and, e.g., the 
Metropolis algorithm with p = 1 is more efficient. 

For the antiferromagnetic initial condition, all spins 
are initially fiippable. The number of fiippable spins de- 
creases rapidly at early times and then the decrease slows 
down as the system coarsens. After each update event, 
the time is incremented by l/# (fiippable spins). Thus 
in one time unit, each fiippable spin changes its state 
once on average. This update step is applied repeatedly 
and the dynamics is averaged over many realizations to 
determine the evolution of the system. 



For times beyond the coarsening time (which scales as 
L^), the energy evolution is characterized by long periods 
where only zero-energy spins can fiip (those with equal 
numbers of up and down neighbors). These long peri- 
ods of stasis are punctuated by progressively more rare 
energy-decreasing events (inset to Fig. [3]). Thus the sys- 
tem typically wanders on successive plateaux that define 
a set of iso-energy points in state space. Occasionally 
there is a drop to a lower plateau by an energy- lowering 
spin-fiip event (Fig. 2]). This feature is illustrated in 
Fig. [5l where we plot the configuration averaged time 
Atn between successive energy-lowering spin-fiip events 
as a function of tn^ the average time at which the n^^ such 
event occurred. Over a substantial range, Atn appears to 
grow exponentially with t^, so that most of the evolution 
is spent wandering aimlessly on iso-energy plateaux. A 
somewhat related continuum picture of this state space 
evolution is presented in Refs. f22l. [23|. 




III. ACCELERATION ALGORITHM 

The evolution of the system becomes so slow that the 
standard Glauber dynamics algorithm described above 
is inadequate to probe the long-time properties of even 
reasonably-sized systems. Figure [3] illustrates this slow- 
ing down for the energy evolution of a typical realiza- 
tion of a 20^ system. The main panel shows the time 
dependence of the gap between the actual energy and 
the ground state energy divided by the total number of 
spins. Henceforth, we term this normalized difference as 
the "energy" El. For this example, the data are roughly 
consistent with the El t'^^^ for 1 < t < 10^. 
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FIG. 3. Time dependence of the energy for a representative 
realization of a 20^ system whose energy stops evolving at 
t ~ 41712. The inset shows the long-time tail of the same 
data; note that the abscissa is Int. 



FIG. 4. Schematic illustration of the state of the system wan- 
dering on fixed-energy plateaux at long times. 




FIG. 5. Average time Atn between successive energy- lowering 
spin-flip events as a function of tn, the time for the n^^ such 
event for 1024 realizations ofalOxlOxlO system. The data 
are smoothed over a 100-point range. 



To reduce the time spent in simulating these iso-energy 
wanderings, we constructed an acceleration algorithm 
and tested that the long-time state achieved by this al- 
gorithm is virtually identical to that of the true zero- 
temperature Glauber dynamics. Our method is based on 
imposing a weak field while the system is on any single 
iso-energy plateau to reduce the time needed to find the 
next energy lowering event. The field is reversed after 
each such event so that the average field is zero. The 
steps of our method are the following: 
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• Apply Glauber dynamics (to fiippable spins only) 
until t = 5L^. This time is sufficiently beyond the 
coarsening time that energy-lowering events have 
become rare (Fig. [3|). 

• For t > 51/^, apply an infinitesimal field to drive the 
state-space motion along a fixed-energy plateau. 
Thus the next energy-lowering event (if it exists) 
is found more rapidly than if no field was applied. 

• After each energy- lowering event occurs, the sign 
of the field is reversed, leading to the alternating 
state-space motion sketched in Fig. HI 

• If the number of active spins goes to zero with- 
out a drop in energy while the field is applied, the 
field is reversed. From this configuration, the sys- 
tem evolves by zero-temperature Glauber dynam- 
ics with the reversed field. If the number of active 
spins again goes to zero without a drop in energy, 
then the final energy value has been reached and 
the simulation is finished. 

We verified that this acceleration algorithm accu- 
rately reproduces the energy obtained by straightforward 
Glauber dynamics for L < 10, where a direct check of this 
acceleration method is computationally feasible. For this 
check, we take all configurations that have been evolved 
to t = by zero-temperature Glauber dynamics and 
evolve each one both by continuing the Glauber dynam- 
ics and by our acceleration algorithm until no fiippable 
spins remain. Over 10^ realizations the fractional differ- 
ence between the energies obtained by these two methods 
is < 1.4 X 10~^. Moreover, the distributions of the final 
energies are virtually indistinguishable. 

For larger L, time to reach the final energy by zero- 
temperature Glauber dynamics is too long to amass suffi- 
cient statistics. Instead, we compare the ultimate energy 
that is reached by starting the acceleration algorithm at 
progressively later times. As shown in the Table |I] in Ap- 
pendix [Al the energy that is ultimately reached changes 
by a negligible amount as the cutoff time is increased. 
For example, for L = 100, the error in the final energy 
reached by the acceleration algorithm is of the order of 
5 X 10~^. Moreover, as illustrated in this table, the ac- 
celeration algorithm is considerably faster than the zero- 
temperature Glauber dynamics. Thus we use the accel- 
eration algorithm for all of our simulation results. On a 
32-core machine, we are able to simulate 10^ realizations 
of a 90^ system in approximately ten days of running 
time. 



IV. THE LONG-TIME STATE 

At sufficiently long times, the energy of each realization 
stops decreasing, either because a blinker configuration 
is reached or occasionally a static final state is reached. 
Figure [6] shows that the probability of reaching a blinker 
configuration approaches 1 as L oo. This observa- 
tion is one of our main results, for which preliminary 



corroboration was given earlier [5, 6j. Somewhat surpris- 
ingly, for the initial condition where the magnetization 
is zero but with otherwise uncorrelated spins, the prob- 
ability to reach a stationary state vanishes more rapidly 
in the L ^ oo limit than for the antiferromagnetic ini- 
tial condition. Even though a blinker consists of a set of 
states of the same energy, we term the "final state" these 
constant-energy configuration (s) that are reached when 
the energy stops decreasing. 
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FIG. 6. Plot of 1 — Pb, the complement of the probability 
Pb of reaching a blinker configuration (o), and 1 — P2, the 
complement of the probability P2 of reaching a state that 
consists of 2 clusters (A), as a function of 1/L. 

In this final state, the spins almost always parti- 
tion into two and only two interpenetrating clusters. 
Using the Hoshen-Kopelman cluster-multilabeling algo- 
rithm to determine the distribution of the number of 
clusters, we find that the probability to find a single clus- 
ter (corresponding to the ground state) or more than two 
clusters in the final state rapidly decays with L (Fig. [6]). 
Final states that contain more than two clusters typi- 
cally consist of multiple narrow filaments of one phase in 
a surrounding background of the opposite phase. In all 
of our simulations on the L x L x L cubes with periodic 
boundary conditions, the largest number of clusters ob- 
served in any realization was seven. (This happened to 
occur in a 38^ system, where the realization consisted of 
six narrow filaments of one phase in a background of the 
opposite phase.) 



A. Relaxation of Blinker Configurations 

As we now discuss, blinker states are responsible for 
the very slow relaxation of the spin system at long times. 
To appreciate the underlying mechanism, it is instruc- 
tive to study the dynamics of the synthetic blinker states 
shown in Fig. [71 In this example, the domain of up 
spins consists of three orthogonal 4 x 12 slabs, each of 
which wraps periodically, so that the apparent slab cor- 
ners are merely visual artifacts. In the cavity defined by 



the confluence of the three slabs, a bhnker state exists. 
By zero-energy spin flips, the portion of the interface 
deflned by this cavity can be in the extremes of fuhy de- 
flated (left panel of Fig. [7|), or fully inflated (right), or 
in some intermediate state (middle). Note that each ex- 
treme conflguration possesses a single blinker spin, while 
intermediate conflgurations have more than one blinker 
spin. Although each blinker spin does not experience any 
energetic bias, there is an effective geometric bias that 
drives the interface to the half-inflated state. This effec- 
tive bias is controlled by the difference in the number of 
flippable spins on the convex and concave corners on the 
interface, N-^ and 7V_, respectively. When the cube is 
mostly inflated — N- is positive so that the interface 
tends to deflate, and vice versa when the cube is mostly 
deflated. Thus the effective bias drives the interface to 
the half-inflated state. 

We quantify the evolution of this synthetic blinker by 
the average flrst-passage time (t) for diii i x i x i fully- 
deflated blinker to become fully inflated. To estimate 
this time, it is helpful to flrst consider the corresponding 
two-dimensional system (Fig. [8]). Near the half-inflated 
state, the interface consists of outer corners and 
inner corners, with N-^ — N- = 1 in two dimensions and 
N-^ ^ i. In a single time unit, all eligible spins on the 
interface flip once, on average. Since — N- = 1, the 
interface area typically decreases by 1. Thus we infer 
an interface velocity u = A A/ At ~ —1. Similarly, the 
mean-square change in the interface area is of the order 
of i y/A. Thus the effective diffusion coeffi- 

cient D is proportional to £. The flrst-passage time is 
dominated by the time to move from the half-inflated 
state to the fully-inflated state by flipping of the order 
oi A = l'^ spins. Since this process is moving against the 
effective bias, the dominant Arrhenius factor ^] in the 
flrst-passage time is r ~ exp(|?i|^^/L)), so that 
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FIG. 8. Two-dimensional analog of the blinker states in Fig. [71 



order of ^; this quantity coincides with the interface ve- 
locity. Similarly, the mean-square change in the interface 
volume in a unit time, which coincides with the diffusiv- 
ity, is of the order of N± ^ D. Consequently, the 

leading behavior of the flrst-passage time is 



Inr - ui^/D 



(4) 



The direct generalization of this argument to d dimen- 
sions gives Inr ~ £^~^. Related aspects of slow evolution 
in three dimensions were discussed for the homogeneous 
kinetic Ising ferromagnet [27] and for a kinetic Ising sys- 
tem with competing ferromagnetic and antiferromagnetic 
interactions [26]. 

Figure [9] shows simulation data for the flrst-passage 
time from the fully-deflated to the fully-inflated state 
in two and three dimensions. The agreement between 
Eq. (|3]) and the two-dimensional data is excellent. In 
three dimensions, simulations are necessarily limited to 
very small ^, while our crude argument is asymptotic; 
nevertheless, the data are qualitatively consistent with 
Eq. dU. The salient result is that the time for a fully- 
deflated blinker to become fully inflated grows rapidly 
with £. 
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FIG. 7. (color online) An 8^ blinker on a 20^ cubic lattice, 
showing the fully-deflated state (left), an intermediate state 
(middle), and the fully-inflated state (right). The bounding 
slabs wrap periodically in the three Cartesian directions. 



For the corresponding three-dimensional system of vol- 
ume V = i^, there are typically N± ~ £^ outer and inner 
corners when the interface is half inflated. The disparity 
in their number is now of the order of i. Thus in a sin- 
gle time step the displacement of the interface is of the 



FIG. 9. First-passage time from the fully-deflated to the fully- 
inflated state for d — 2 (o) and d — 3 (A). The line for d — 2 
is the best fit r = 1.40 exp(1.33^), while the curve for d = 3, 
r = 3 exp(0.8^^), is merely a guide for the eye. Each data 
point is based on at least 128 realizations. 



From the dynamics of the synthetic blinker of Fig. [71 we 
can now understand the long time scales associated with 
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the relaxation of a large system. Indeed, consider two 
such blinkers that are oppositely oriented and spatially 
separated so that they do not overlap when both are 
half inflated, but just touch corner to corner when both 
are inflated (Fig. [TQ|) . As long as the blinkers do not 
overlap, their fluctuations do not change the energy of 
the system. However, if these blinkers touch, then a spin 
flip event has occurred that lowers the energy. After this 
irreversible coalescence, subsequent spin flip events cause 
the two blinkers to ultimately merge. 





FIG. 10. A schematic two-dimensional projection of a blinker 
coalescence event. The tips of the blinkers correspond to the 
fully-inflated configuration shown in the right panel of Fig. [71 



spin. This energy decreases with L in a manner consis- 
tent with the power-law dependence El ~ (Fig. [TT]). 
However, there is systematic curvature in this data, and 
we extrapolate the local two-point slopes in the plot of 
El versus L to obtain the estimate e ^ 1, in agreement 
with previous results based on smaller-scale simulations 
[6] . This dependence implies that the total interface area 
between spin domains scales as L^. 
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Because of the impermanence of a two-blinker config- 
uration, we term it a pseudo-blinker. We assert that 
each coalescence of a pseudo-blinker corresponds to one 
of the increasingly rare energy-lowering spin-flip events 
sketched in Fig. |4j The time for pseudo-blinker coales- 
cence is extraordinarily long because the time for each 
blinker to reach a nearly-inflated state is a rapidly in- 
creasing function of its size. These coalescences corre- 
spond to energy- lowering spin- flip events at long times. 

We now turn to the true blinker states. As might be 
anticipated from the example in Fig. [TJ simulations indi- 
cate that the fraction of blinker spins is small — of the 
order 3 x 10"^ to 4 x 10"^ of all spins for systems of linear 
dimension L < 50. The instantaneous number of blinker 
spins also fluctuates substantially so that their number is 
not a meaningful characteristic of the blinker states. A 
more robust measure is the total volume that is accessed 
by blinker spins. We determine this accessible volume as 
follows: Once a true blinker state is first reached (which 
we define as !Bo), we drive the system with an infinites- 
imal positive field until the spin configuration with 
no flippable spins, is reached. Then starting again from 
^0, we drive the system with an infinitesimal negative 
field until there are no flippable spins and the configu- 
ration ^_ is reached. The difference — defines 
the total blinker volume. The resulting blinker volume 
fraction is a slowly increasing function of L and extrap- 
olating to L ^ oo gives an asymptotic blinker volume 
fraction of approximately 9% — a finite fraction of the 
entire system. 

B. Asymptotic Energy 

An important characteristic of a finite system of lin- 
ear dimension L is its energy El infinite time. As 
mentioned in Sect. IIIIl what we term the energy is ac- 
tually the energy gap above the ground state energy per 



FIG. 11. Normalized average energy (o) and genus of the final 
state (A) as a function of L. The relative error for each data 
point is less than 1.4 x 10~^. 
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FIG. 12. The normalized final- state energy distribution for 
L = 54 (o), L = 76, (A), and L = 90 (y). 



The scaled distribution of energies, P{E/El)^ exhibits 
an excellent data collapse (Fig. [12]). The distribution has 
a well-defined peak that is close to a Gaussian; there is 
also a noticeable linear tail at low energies. The fact that 
the energy distribution P{E/El) remains broad in the 
thermodynamic limit is not surprising as all our numeri- 
cal findings show the lack of self- averaging. It would be 
interesting to understand qualitatively the shape of the 
energy distribution P{E/El)^ or at least the asymptotic 
behaviors in the E/El ^ and E/El oo limits. 
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C. Domain Topology 

The long-time state of the three-dimensional system 
is, in general, topologically complex because it is possi- 
ble — in fact likely — that the interface contains many 
holes [28|. The number of holes is also known as the 
genus g. Formally, the genus of a connected, orientable 
surface is an integer that represents the maximum num- 
ber of cuts that can be made through the surface along 
non-intersecting closed simple curves without disconnect- 
ing the resulting manifold. As elementary examples, the 
genus of a sphere is ^ = 0, while the genus of a doughnut 
is ^ = 1. 

FIG. 13. (color online) Simple examples of interfaces with 
genus g — 2 and g = S for a periodic 2x2x2 system. 




FIG. 14. (color online) Examples of a low-genus (g = 3) and 
high-genus (g = 16) domain on a 20^ lattice. 

The genus of a surface can be expressed through its 
Euler characteristic x (see e.g. (29l, [3Q|) 

X = 2{l-g) = V-E + J. (5) 

The latter equality relates the Euler characteristic to 
easily-measured features of the domain interface: V, 
the number of vertices on the interface, £, the number 
of edges, and S^, the number of faces. As simple ex- 
amples, the Euler characteristic of an isolated cube is 
X = 8 — 12 + 6 = 2, corresponding to genus 0. (Topolog- 
ically, the cube is identical to the ball, so the boundary 
of the cube is a sphere.) The Euler characteristic of a 
linear filament that wraps around the torus in one di- 
rection is zero. By discretizing this filament as a 2 x 1 
cluster that wraps onto itself, we have V = 8, £ = 16 and 
9^ = 8, corresponding to genus g = I. Note that the Eu- 
ler characteristic does not depend on the length scale of 
the discretization. Similarly for a cluster that percolates 
in two directions (Fig. [T3l) . y = 8 - 20 + 10 = -2 so the 



genus g = 2. Finally for a cluster that percolates in all 
three Cartesian directions, x = 8 — 24 + 12 = —4 so the 
genus g = 3. 

Since blinker spins do not affect the topology, we freeze 
these spins in their orientations at the time when the in- 
terface topology is measured. To measure the topology, 
we first identify all the clusters in the final state by the 
Hoshen-Kopelman algorithm [24]. We then compute the 
Euler characteristic of each cluster. If there are only two 
clusters, then by construction they have the same inter- 
face and thus the same Euler characteristic. If a final 
state has more than two clusters, we use the maximum 
genus among all clusters as the genus of the system. To 
determine the Euler characteristic numerically, we first 
determine the number of faces 3^. This quantity equals 
the number of neighboring antiparallel spins (accounting 
for the periodic boundary conditions) which, in turn, is 
directly related to the energy of the system. Once we 
identify a new face on the interface, each of the four ver- 
tices and the four edges bounding this face are added to 
the current counts of V and £, as long as they have not 
yet already been counted as part of another previously- 
encountered face. 




g/<g> 



FIG. 15. The final-state genus distribution for L = 54 (o), 
L = 76, (A), andL = 90 (y)- 

The resulting domain topologies are quite diverse 
(Fig. [14]). For example, for L = 20, the smallest genus 
observed was 1, the largest genus was 18, while the av- 
erage genus is approximately 4.36. The average genus 
for a given L, defined as ^l, again appears to grow as 
a power law in L, gL ^ but with substantial finite- 
size corrections (Fig. [11]). Analogous to the behavior for 
the average energy, the data for gL versus L on a double 
logarithmic scale are systematically curved upward and 
extrapolating the effective exponent to L ^ oo yields the 
estimate 7 ~ 1.7. The scaled genus distributions at long 
times for L = 54, 76, and 90 also show excellent data 
collapse (Fig. [T5]) . 

Although we separately studied the energy and the 
genus, these two quantities can be intimately connected 
by simple topological considerations. We start by sim- 
plifying Eq. dSJ. For a closed surface that defines the 
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FIG. 16. A portion of an interface to illustrate vertices (la- 
belled) that are shared among 3, 4, 5, or 6 edges. 

interface between spin domains, each face has 4 edges 
and each edge must be shared between 2 adjacent faces. 
Hence 

£ = 29^. (6) 

Similarly, each edge has 2 vertices and each vertex is 
shared among 3, 4, 5, or 6 adjacent edges (Fig. [T6|). This 
leads to the bounds 

F 2F 

- < V < — . (7) 

3 - - 3 ^ ^ 

Using these two relations in Eq. (j5j) gives 

3" 3" 
-3<X<3^ (8) 

While the lower bound is useful, the upper bound can be 
replaced by the much stronger condition x ^ 2, with the 
maximal value of 2 being achieved for a sphere. With 
this replacement, we obtain the following bounds on the 
genus: 

0<^<| + 1. (9) 

However, 1 is directly related to the total energy of the 
system because each face corresponds to a single pair of 
anti-aligned spins per spin and thus has an energy cost of 
+2 (when the interaction strength is set to 1, see Eq. ([T])). 
Thus 3" = L^El. 

Finally, we make the guess that the actual value of 
g lies roughly midway between the upper and lower 
bounds; namely, g oc 3^. Using our exponent definitions 
El ~ and (g) ~ L^, our argument leads to the ex- 
ponent inequality 

e + 7<3. (10) 

Our numerical estimates for these two exponents given 
above, e ~ 1 and 7 ^ 1.7 are consistent with Eq. (p!Q|) . 

We can carry this analysis a bit further by making 
use of some simple facts in discrete differential geome- 
try 0. Let us define the number of vertices with m 
incident edges as V^. It is useful to introduce the notion 
of a "defect" that is associated with each vertex. The 
defect for a vertex is defined as the difference between 



the sum of the angles of all the faces at the vertex and 
27r. With this definition, it is easy to see that the defect 
of a vertices of types 3, 4, 5, and 6 are ^, 0, — and 
— TT, respectively. For any domain, the total defect of all 
vertices on the surface equals 27rx; this is essentially the 
Gauss-Bonnet theorem for a discrete interface j33| . Thus 
we have the general relation (see Eq. (|5])) 

I V3 - I V5 - ^ Ve = 2^ X = 4^(1 - (11) 

Therefore the genus of a surface and the number of ver- 
tices of various types are related by 

5=l + i(2V6 + V5-V3). (12) 

From the examples of domains shown throughout this 
work, almost all vertices are of type m = 4 while vertices 
of type m = 5 seem to be next most common. Vertices 
of type m = 3 are associated with blinkers and therefore 
should be few in number. Vertices of type m = 6 arise at 
a 3-fold branch of the domain and therefore seem to be 
the most rare. Thus we expect that vertices of type m = 
5 scale the same way as (g) which numerically appears 
to grow as L^-^ . 

D. Survival Probability 

The relaxation process is naturally characterized by 
S{t)^ the probability that the energy of the system is still 
decreasing at time t. Since energy- lowering spin flips oc- 
cur rarely at long times, it is not immediately evident 
whether the most recent energy-lowering spin flip is the 
last such event or whether another energy-lowering flip 
event will occur sometime in the distant future. To de- 
termine if the energy has reached its final value in an 
efficient way, we use the following algorithm that is a 
variant of our acceleration algorithm. As a preliminary, 
we separately track both positive-energy and zero-energy 
flippable spins; the former are those for which the en- 
ergy decreases if such a spin actually flips. We start the 
simulation by running zero-temperature Glauber dynam- 
ics until no positive-energy flippable spins remain. At 
this time, defined as Tq, the configuration Co may have 
reached the final value of the energy. 

We now proceed as follows: 

• Apply an infinitesimal magnetic field. If an energy- 
lowering spin-flip occurs, then the energy of Co is 
not the final energy. In this case, the system is re- 
turned to Co and subsequently evolves by Glauber 
dynamics until again no positive-energy spins re- 
main and a new candidate final configuration and 
final time is reached. 

• If an energy- lowering spin- flip does not occur, the 
system is returned to Co and a field is applied in 
the opposite direction. Again, if an energy- lowering 
spin-flip occurs, the system is returned to Co and 
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10° 10' 10' 10' ^ 10^ 



FIG. 17. Survival probability S{t) versus time t for L = 
4, 6, 8, 10, 14, 20, 30 and 40 (lower left to upper right) on a 
double logarithmic scale. The averaging has been performed 
over 10^ realizations for L < 10, 10240 realizations for L = 
14, 20, 30, and 2048 realizations for L = 40. 




FIG. 18. S{t) versus Int for L = 20,30, and 40 on a double 
logarithmic scale. 

subsequently evolves by Glauber dynamics until no 
positive-energy spins remain and a new candidate 
final state is reached. 
• If an energy-lowering spin-flip does not occur after 
the field has been applied in both directions, then 
Co is at the final energy and the survival time equals 
To- 

We record the time when the energy of a system stops 
changing, from which we infer the time dependence of 
the survival probability S{t). This time dependence is 
both surprisingly complex and extremely slow for small 
system sizes (Fig. [T7|). For example, for L = 10, 40 re- 
alizations had not yet relaxed to their ultimate energy 
by t = 10^, a time that is seven orders of magnitude be- 
yond the coarsening time scale of 10^. The data cutoff 
for L > 10 was imposed because of CPU time limitations. 
By L = 20, the dependence of S{t) on t becomes smooth 
and reasonably systematic and a plot of S{t) versus Int 
on double logarithmic scale (Fig. [T8|) suggests that the 
long-time data can be reasonably fit to an inverse loga- 



rithmic dependence S{t) ~ (Int) ^, with a ^ 3. 

V. DISCUSSION 

We investigated the evolution of the kinetic Ising model 
that is endowed with single spin-flip dynamics on a finite 
cubic lattice with periodic boundary conditions. The sys- 
tem starts in the antiferromagnetic state and is quenched 
to zero temperature. The details of the initial condi- 
tions are secondary as long as the magnetization van- 
ishes. (If the initial magnetization is non-zero, the evolu- 
tion is much simpler and the Ising ferromagnet falls into a 
ground state.) We asked the simple question: what hap- 
pens? A natural expectation might be that the ground 
state should be reached. A more comprehensive version 
of this presumption is encapsulated by the central dogma 
of coarsening, which asserts: 

1. Ising ferromagnets have just two metastable states, 
which coincide with the ground states. 

2. If an Ising ferromagnet is endowed with zero- 
temperature single spin- flip dynamics, or more gen- 
erally with a non order-parameter conserving dy- 
namics, then one of the two ground states is neces- 
sarily reached. 

3. The time to reach a ground state scales with the 
linear dimension of the system as L^. 

This central dogma is indeed correct in one dimension. 
However, for two-dimensional Ising ferromagnets, there 
are numerous metastable states that consist of single- 
phase stripes whose total number grows as M2 ~ Q, 
where g = ^(V^+1) is the golden ratio. Nevertheless, the 
failure of the central dogma in two dimensions is rather 
benign, as one of the ground states is reached [5.,£.,^] with 
probability close to 2/3. Moreover, for most realizations, 
the final state (either a ground state or a stripe state) is 
approached in a time that scales as L^. 

In three dimensions, however, the central dogma com- 
pletely fails; viz., all its three basic tenets are wrong. 
First, the number of metastable states M3 scales expo- 
nentially with the system size: InMs ~ [31]. Further, 
the ground states are never reached (for sufficiently large 
systems) and the relaxation time is anomalously long. 
We provided heuristic and numerical evidence that the 
relaxation time scales as . Thus for a macroscopic 
system with L ~ 10^ the relaxation time considerably 
exceeds any time scale in the Universe. 

Since the approach to the long-time state is extraordi- 
narily slow, even for a system as small as 10 x 10 x 10, 
there are still realizations (albeit a small fraction) for 
which the energy has not yet reached its final value by 
t = 10^, whereas the standard coarsening time is of the 
order of 10^. We constructed a physical picture, based 
on the coalescence of blinker-like configurations, that in- 
stead predicts a relaxational time scale that grows ex- 
ponentially in the linear dimension of the system. In 
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particular, the survival probability S'(t), defined as the 
probability that the energy has not yet relaxed to its fi- 
nal value by time t seems to decay as a power law in 
1/lnt. While the mechanism of blinker coalescence ap- 
pears plausible, we do not have a theoretical explanation 
for the functional form of S{t). The primary feature of 
the relaxation that we wish to emphasize is that the stan- 
dard picture of coarsening, characterized by a time that 
scales as L^, is inappropriate for the three-dimensional 
Ising model with zero-temperature Glauber dynamics. 

Another striking feature of the three-dimensional Ising 
ferromagnet is that a set of connected metastable mi- 
crostates with equal energy is reached rather than a 
'frozen' metastable state. Each such set of states con- 
tains a small number number of flippable blinker spins 
that can flip ad infinitum without any energy cost. Thus 
the system can wander within one of these iso-energy 
sets forever. Even though the number of blinker spins 
is small, the spatial volume over which the blinker spins 
can roam comprises of the order of 10% of the volume of 
system in the limit of large L. 

The topology of the long-time state is much richer than 
that of the corresponding two-dimensional system. In 
two dimensions, the only possible states at infinite time 
are the ground state or an even number of alternating 
single-phase stripes (for periodic boundary conditions). 
In contrast, the long-time states in three dimensions are 
highly interpenetrating and contain many holes fFig.[T4|). 
Correspondingly, the average genus of the domains scales 
with linear dimension as L^-^. Aside from this global 
characterization of the domains, it is not clear what are 
the most useful measures of the domain geometry. 

While we focused on the cubic lattice, similar behav- 
iors should arise for the kinetic Ising model on other 
e^;en-coordinated lattices in three dimensions. We delib- 
erately avoided ot/tZ-coordinated lattices or other complex 
networks where the coordination can be odd, as the zero- 
temperature Ising- Glauber system quickly freezes, in dis- 
agreement with the central dogma predictions. However, 
this freezing has a local and trivial nature. For example, 
for an odd-coordinated lattice, single-phase droplets can 
arise, in which spins within a droplet each have more in- 
ternal than external neighbors [32j, l34H37i] , as illustrated 
for the example of a hexagonal droplet on the hexag- 
onal lattice (Fig. [T9|) . These droplets can thus remain 
forever in the phase opposite to that of the background. 
Another example of this local freezing is the Ising model 
with zero-temperature Kawasaki (spin exchange) dynam- 
ics where again local defects quickly arise that stop 
the overall relaxation process. 

Intriguing and mostly unexplored behaviors arise for 
non-cubic systems (Fig. [20]), for example, a L x L x aL 
system. When the aspect ratio a is small, the system be- 
comes a thin square slab. We are generically interested 
in the thermodynamic limit L ^ oo with a fixed, so a 
thin square slab does not reduce to a two-dimensional 
system. When the slab is thin, the long-time state re- 
sembles Swiss cheese, with directed holes perpendicular 




FIG. 19. Example of a frozen cluster of up spins on the hexag- 
onal lattice. 

to the slab. Because of the periodic boundary condi- 
tion in all directions, there is no possibility of forming 
the stripe states that arise in the two-dimensional sys- 
tem p\ia\. For a 1, corresponding to a long bar, the 
long-time state consists of a series of alternating domains 
of the two phases. As the bar becomes wider, percolation 
in the long direction eventually occurs, and the geome- 
try begins to resemble the plumber's nightmare (middle 
panel of Fig. [2Q|) . 




FIG. 20. Example long-time states for a 32 x 32 x 8 slab, a 
32^ cube, and a 8 X 8 X 32 bar. 

It is striking that the three-dimensional system is so 
much more complicated than the corresponding two- 
dimensional case. This fact suggests that many more 
surprises await discovery for the kinetic Ising model in 
higher dimensions. Moreover, most of our findings are 
empirical in nature and they beg for the development of 
new theoretical perspective and geometrical descriptions 
of the domains. Another challenging extension of the 
present work is to describe the fate of spin systems with 
a non-scalar order parameter, such as the XY model, that 
are quenched to zero temperature 

We gratefully acknowledge financial support from NSF 
grant DMR0906504 (JO and SR) and NSF grant GGF- 
0829541 (PLK). 



Appendix A: Accuracy of the Acceleration 
Algorithm 

To test the accuracy of our acceleration algorithm, 
the system is evolved until a cutoff time r with zero- 
temperature Glauber dynamics. Then at time r, an in- 
finitesimal magnetic field is applied, that alternates as 
the system descends successive energy plateaux, until no 
flippable spins remain. Table Ugives the final energy that 
is reached (where no flippable spins remain) when the ac- 
celeration algorithm is applied for different cutoff times 



r. These data are shown for the cases L = 10, 20, and 
100, with 10^, 10240, and 128 reahzations, respectively. 
The extremely weak dependence of the final energy as a 
function of r indicates the level of accuracy of the accel- 
eration algorithm. 



L 


7" 


(Er) 


R 


10 


500 


.4245900960 






10^ 


.4245901020 


34 




10^° 


.4245901020 


60 


20 


2000 


.283285352 






10^ 


.283287939 


1.4 




10^ 


.283287939 


5.2 




10^ 


.283287939 


39 




10^ 


.283288232 


378 


100 


5 X 10^ 


.083666469 






10^ 


.083663406 


1.2 




5 X 10^ 


.083662656 


4.7 




10^ 


.083662594 


7.2 




10^ 


.083662531 


67 



TABLE I. Average final energies {El) for different cutoff 
times T and system sizes L. 



The last column gives the ratio R of the CPU time 
needed to simulate the system until no fiippable spins re- 
main when the acceleration algorithm is imposed at the 
cutoff time r compared to imposing this algorithm at 
r = hi? . For example, it took 67 times longer to sim- 
ulate the L = 100 system by running zero-temperature 
Glauber dynamics to t = 10*" and subsequently impos- 
ing the acceleration algorithm compared to running zero- 
temperature Glauber dynamics to t = 5 x 10^ and then 
imposing the acceleration algorithm. The relative differ- 
ence in the energies by the two protocols is approximately 
5 X 10~^, thus providing justification for our use of the 
acceleration algorithm at t = hi? . 

Appendix B: Small Systems, Blinker States, 
Number of Clusters 

The evolution of the smallest possible lattice L = 2 
helps to illustrate the complexities of larger systems. 
When L = 2, there are 8 spins and 2^ = 256 possible 
states that can be enumerated to determine all details 
of the evolution. There are also only two possible final 
states: the ferromagnetic ground state (F) and a static 
metastable state (M) that consists of a square of four 
spins of one sign and an adjacent square of four spins of 
the opposite sign. There are nine distinct paths in state 
space that start at the antiferromagnetic state and end 
at these two final states (Table The average survival 
time until the system reaches one of the final states is 
11^ = 1.841666..., while the probability of ultimately 
reaching the F final state is 
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path 


# flips 


time 


prob. 


final state 


1 


4 


43 
56 


2 

21 


M 


2 


4 


143 
168 


1 

14 


M 


3 


6 


341 
280 


1 

21 


M 


4 


4 


85 
56 


3 

14 


F 


5 


6 


1583 
840 


1 

35 


F 


6 


6 


1793 
840 


4 

35 


F 


7 


6 


127 
56 


4 
21 


F 


8 


6 


395 
168 


1 
7 


F 


9 


8 


761 


2 


F 




280 


21 





TABLE II. The nine state-space paths to the final state for 
L — 2 starting from the antiferromagnetic initial condition. 
Listed for each path are the number of spin flips until the flnal 
state is reached, the time to reach the flnal state on the path, 
and the probability of the path. Also listed is the flnal state 
for each path, either metastable (M) or ferromagnetic (F). 



A complete enumeration is already not feasible for lin- 
ear dimension L = 4, where the number of states is 
2^^ ^ 1.84 X 10^^; thus simulations are necessary when 
L > 4. For L = 4, the average survival time is now 6.16, 
while the longest survival time observed in any realization 
is 22.2. For L = 6 and 8, the respective average survival 
times are 11.9 and 27.9, while the longest survival times 
are 128 and 3.43 x 10^. For L = 10, realizations that live 
beyond t = 10^^ are possible, although rare, and it is not 
possible to quote an average survival time. The existence 
of such long-lived realizations for L > 10 contributes to 
the difficulty in the understanding of the behavior of the 
survival probability. 



L 




Pf 


Pb 


2 


11 

14 


3 

14 





4 


0.6814(1) 


0.3186(1) 





6 


0.3523(2) 


0.6353(2) 


0.01246(4) 


8 


0.1842(1) 


0.7373(1) 


0.07853(9) 


10 


0.1045(1) 


0.7170(1) 


0.1785(1) 


20 


0.01377(4) 


0.3059(1) 


0.6803(1) 


32 


0.00322(8) 


0.1091(4) 


0.8877(4) 


54 


0.00066(6) 


0.0406(4) 


0.9587(4) 


76 


0.00040(6) 


0.0250(5) 


0.9746(5) 


90 


0.00039(6) 


0.0199(4) 


0.9797(4) 



TABLE III. Probabilities of reaching the ground state Pg, the 
frozen state P/, and the blinker state Pb versus L. The error 
in the last digit is shown in parentheses. 



It is also worth noting that infinitely long-lived blinker 
states first appear for the case of L = 5. Table IIIII 
gives the data for the probabilities of reaching the ground 
state, a frozen static state, or a blinker state as a func- 
tion of L. The former two probabilities decrease rapidly 
with L and appear to approach zero for large L, while the 
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probability that a blinker state is reached approaches 1 
as L increases (see also Fig. [6]). As mentioned at the 
beginning of Sect. |lVl the long-time state almost always 
consists of exactly two clusters. Table HVl gives the proba- 
bilities that the long-time state consists of one, two, three 
or more than three clusters as a function of L. 



L 


P(l) 


P(2) 


P(3) 


P(>3) 


2 
4 
6 
8 

10 
20 

32 
54 
76 
90 


11 

14 

.6814(1) 
.3523(2) 
.1842(1) 
.1045(1) 
.01377(4) 
.00322(8) 
.00066(6) 
.00040(6) 
.00039(6) 


3 

14 

.3186(1) 
.6475(2) 
.8128(1) 
.8866(1) 
.96052(6) 
.9720(2) 
.9802(3) 
.9824(4) 
.9839(4) 





.000245(5) 
.00303(2) 
.00893(3) 
.02475(5) 
.0230(2) 
.0171(3) 
.0150(4) 
.0133(4) 




<. 0000001 
<. 0000001 
.0000004(2) 
.000015(1) 
.00096(1) 
.00180(6) 
.0020(1) 
.0022(1) 
.0024(2) 



TABLE IV. Probabilities of reaching a state that contains 1, 
2, 3, or > 3 clusters at long times, with the error on the last 
digit in parentheses. 
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